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■  Abstract  We  review  the  use  of  ray  models  for  internal  waves,  particularly  formu¬ 
lations  for  calculating  wave  amplitudes  along  the  ray.  These  are  expressed  in  spatial, 
wave  number,  and  phase-space  coordinates.  The  choice  of  formulation  affects  not  only 
the  difficulty  of  the  calculations  for  rays  and  caustics  but  also  the  degree  to  which 
the  waves  satisfy  slowly  varying  assumptions.  We  describe  several  examples  taken 
from  atmospheric  and  oceanic  applications  that  illustrate  the  variety  of  options  for  ray 
models. 


1.  INTRODUCTION 

Ray  models  are  the  basis  for  much  of  our  understanding  of  atmospheric  and  oceanic 
internal  waves.  They  describe  the  wave  field  that  emerges  from  various  sources, 
the  subsequent  propagation  through  a  nonuniform  background,  and  the  approach 
to  dissipation.  Ray  models  also  represent  some  of  the  leading  attempts  to  account 
for  internal-wave  spectra  and  internal-wave  dissipation  rates  in  parts  of  the  atmo¬ 
sphere  and  ocean,  and  they  have  been  used  to  parameterize  internal-wave  drag 
in  atmospheric  circulation  models.  Here  we  review  and  relate  ray  models  of  in¬ 
ternal  waves  for  these  and  other  applications.  We  concentrate  on  the  formulation 
rather  than  on  their  predictions,  a  subject  included  in  other  surveys  (e.g.,  Fritts  & 
Alexander  2003). 

Ray  models  can  be  formulated  in  spatial  coordinates,  in  wave-number  coor¬ 
dinates,  or  in  a  mix  of  the  two.  The  formulation  in  wave-number  coordinates  is 
equivalent  to  the  ray  description  of  the  Fourier  transform  of  the  wave  field.  How¬ 
ever  formulated,  the  wave  amplitudes  are  controlled  by  the  same  basic  rules  of 
ray  convergence  and  divergence,  but  some  formulations  are  more  convenient  than 
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others.  For  example,  for  a  height-dependent  horizontally  uniform  background,  the 
transformation  from  horizontal  spatial  coordinates  to  horizontal  wave-number  co¬ 
ordinates  straightens  the  rays  and  eliminates  the  more  complicated  occurrences  of 
the  ray  singularity  known  as  a  caustic,  so  the  ray  calculation  and  the  correction  of 
the  caustic  can  be  made  with  relative  ease. 

A  caustic  is  a  singularity  particular  to  ray  theory.  It  occurs  where  neighbor¬ 
ing  rays  intersect  each  other,  resulting  in  an  extreme  breakdown  of  ray  theory’s 
slowly  varying  approximation  and  a  ray  prediction  of  infinite  wave  amplitudes. 
The  familiar  case  is  the  buoyancy-frequency  turning  point  (Lighthill  1978,  section 
4.11),  where  neighboring  rays  intersect  as  they  reverse  their  vertical  direction  of 
propagation.  This  case  is  for  a  particular  model  and  for  a  particular  formulation 
of  the  ray  equations.  Typically,  the  caustic  locations  are  more  widespread,  and  in 
some  models  they  can  occur  almost  anywhere  along  the  ray. 

The  problem  with  caustics  is  that  their  correction  in  numerical  ray  tracing  is 
generally  nontrivial.  Very  few  of  the  models  that  we  review  correct  any  caustics  at 
all.  For  some  applications,  such  as  those  concerned  with  the  approach  to  a  critical 
layer,  caustics  are  of  limited  interest.  The  rays  may  pass  through  a  caustic  on  the 
way  to  a  critical  layer,  but  assuming  no  dissipation  takes  place  at  the  caustic,  it  is 
enough  to  know  that  the  amount  of  wave  action  carried  by  the  waves  is  conserved 
through  the  caustic,  and  that  ray  theory  becomes  valid  again  after  the  ray  leaves 
the  caustic.  In  applications  where  caustics  occur  at  locations  of  interest,  the  choice 
of  model  formulation  can  ease  or  obviate  the  correction  of  caustics.  Part  of  the  aim 
of  this  paper  is  to  review  how  this  has  been  done. 

We  start  in  Section  2  with  a  description  of  the  spatial  and  wave-number  formu¬ 
lations  of  ray  theory.  Initially  we  assume  a  steady  source  of  waves  in  a  steady  back¬ 
ground,  as  do  many  ray  models  of  internal  waves,  but  the  general  time-dependent 
ray  equations  for  wave-amplitude  calculations  are  also  of  interest  and  are  described 
later  in  the  section. 

In  Section  3  we  discuss  caustics.  The  familiar  caustic  at  a  buoyancy-frequency 
turning  point  is  one  case  that  can  be  handled  easily,  by  a  matching  method  or 
by  a  uniform  approximation,  usually  involving  an  Airy  function.  The  details  are 
given  in  many  other  references  (e.g.,  Kravtsov  &  Orlov  1999,  Lighthill  1978),  so 
we  concentrate  on  another  interesting  example  of  a  caustic  that  is  less  familiar  but 
important  in  certain  models  of  internal  waves  generated  by  flow  over  topography  or 
by  an  oscillating  source.  Miles  (1969)  and  Lighthill  (1978)  suggested  alternatives 
to  the  basic  ray  method  that  avoid  the  difficulties  in  these  respective  examples. 
Miles’s  alternative  is  related  to  Maslov’s  method,  the  subject  of  Section  4. 

Maslov’s  method  takes  advantage  of  the  fact  that  the  occurrence  of  caustics  is 
formulation  dependent.  A  caustic  cannot  occur  along  the  ray  in  its  spatial  formu¬ 
lation  and  in  its  wave-number  formulation  at  the  same  time.  The  idea  in  Maslov’s 
method  is  to  calculate  the  ray  solution  in  the  formulation  without  the  caustic,  and 
then  to  map  it  to  a  solution  in  the  other  formulation  by  Fourier  transform.  In  theory, 
this  corrects  all  types  of  caustics,  and  in  some  cases  Maslov’s  method  is  actually 
easy  to  implement.  We  give  an  example  in  Section  4. 
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In  Section  5  we  discuss  a  range  of  applications  that  indicates  the  variety  of  op¬ 
tions  used  in  ray  models  for  internal  waves.  Sometimes  the  formulation  is  chosen 
for  the  variables  of  direct  interest,  but  the  choice  is  also  influenced  by  computa¬ 
tional  limitations  and  concerns  about  caustics.  We  do  not  cover  basic  internal-wave 
theory,  which  can  be  found  in  texts  such  as  Gossard  &  Hooke  (1975),  Lighthill 
(1978),  Gill  (1982),  and  Nappo  (2002).  For  other  reviews  of  internal  waves,  see 
Fritts  &  Alexander  (2003),  McIntyre  (2001),  Baines  (1995),  Wurtele  et  al.  (1996), 
and  Muller  et  al.  (1986). 


2.GENERALTHE0RY 

Lighthill  (1978,  section  4.6)  discusses  the  ray  tracing  of  internal  waves  through 
a  horizontal,  vertically  varying  wind  U(z ),  directed  along  the  x—  axis.  A  ray  is 
defined  as  the  position  x(f)  =  (x(t),  y(t),  z(t ))  that  moves  through  the  medium 
at  the  local  group  velocity  of  the  waves  cg  =  (cgi,  cg 2,  cg 3)  measured  by  a  sta¬ 
tionary  observer  on  the  ground.  Because  U  depends  only  on  z,  the  horizontal 
wave  numbers  k,  l  remain  constant  following  the  ray,  as  does  the  frequency  co 
measured  by  the  stationary  ground  observer.  The  vertical  wave  number  m  varies 
along  the  ray,  as  does  the  intrinsic  frequency  So  =  co  —  kU  measured  by  an 
observer  moving  with  the  local  wind  velocity.  For  a  single  wave  train  of  fixed 
k,  l,  co  and  height-dependent  m(z),  the  wave  amplitudes  are  predicted  from  the 
constancy  of  the  wave-action  flux  Cg^E/cb,  where  E  is  the  wave-energy  den¬ 
sity  (see  Equation  19).  This  simple  case  explains  aspects  of  phenomena  such 
as  the  Booker-Bretherton  critical-layer  interaction  (Gossard  &  Hooke  1975).  To 
review  other  ray  models  of  internal  waves,  we  need  to  consider  the  more  general 
theory. 

The  dispersion  relation  is  &>(x,  t)  —  G(k,  x,  /),  for  time  t,  position  x  =  (x,  y, 
7),  and  wave-number  vector  k  =  (k,  /,  m).  The  wave  frequency  is  co  —  So  +  k  •  U, 
for  intrinsic  frequency  So  and  a  background  velocity  U,  which  can  vary  in  space 
and  time.  It  is  convenient  to  follow  Hayes’s  notation  (1970)  for  specifying  the 
independent  variables.  When  k  and  x  are  both  treated  as  independent  variables, 
we  use  G  and  partial  derivatives  denoted  by  subscripts.  When  k  and  x  are  treated 
as  functions  of  one  another,  we  use  00  and  partial  derivatives  denoted  by  3/3 f, 
V  =  (d/dx,  3/3 y,  3/3 z),  and  Vk  =  (3/3 k,  3/3/,  3/3 m).  Notation  such  as  Gkk 
refers  to  the  tensor  with  components  G,y.,  G/:/,  etc. 

The  ray  equations  are  (Lighthill  1978,  Hayes  1970) 


dx/dt  =  Gk 

(i) 

X 

Cl 

1 

II 

^3 

(2) 

where  d/dt  =  3/3f  +  c?  ■  V  and  c,,  —  Gk  is  the  group  velocity  vector.  There 
is  no  need  yet  to  specify  the  form  for  G,  though  for  now  we  assume  that  G  is 
independent  of  time.  We  relax  this  assumption  later. 
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The  ray  equations  for  wave-amplitude  calculations  can  be  formulated  in  spatial 
coordinates  or  in  wave-number  coordinates.  We  first  describe  how  these  formu¬ 
lations  are  related  to  the  ray  phase-space  formulation,  where  the  independent 
coordinates  are  both  x  and  k. 

In  phase  space,  the  initial  conditions  uniquely  determine  each  ray  path,  so  there 
are  no  ray  intersections  and  hence  no  caustics.  Furthermore,  the  density  of  wave 
action  in  phase  space  is  constant  along  the  ray  (e.g.,  Hertzog  et  al.  2002).  Although 
the  absence  of  caustics  and  the  constancy  of  the  wave-action  density  along  the  ray 
are  nice  simplifications,  we  are  not  usually  interested  in  the  phase-space  density 
of  wave  action  itself  but  rather  in  the  spatial  density  or  wave-number  density  of 
wave  action.  To  obtain  these  densities,  project  the  ray  solution  from  phase  space 
to  the  spatial  domain  or  to  the  wave-number  domain.  In  the  projected  domains, 
two  neighboring  rays  can  project  onto  the  same  point,  resulting  in  a  caustic. 

Figure  1  illustrates  the  case  where  the  projected  rays  form  a  caustic  in  the  spatial 
domain.  It  is  also  possible  to  have  a  caustic  in  the  wave -number  domain,  but  an 


Figure  1  Rays  in  phase  space  are  shown  in  this  schematic  diagram,  along  with  their 
projections  onto  the  spatial  x  domain  and  the  wave-number  k  domain.  Rays  never 
intersect  in  phase  space.  Ray  intersections,  and  hence  caustics,  occur  only  in  the  ray 
projections,  in  this  case  in  the  x  domain. 
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important  point  is  that  a  ray  in  phase  space  cannot  simultaneously  project  onto  a 
caustic  in  the  wave -number  domain  and  onto  a  caustic  in  the  spatial  domain.  At 
a  caustic  in  the  spatial  domain,  neighboring  rays  have  the  same  x  but  different  k. 
The  situation  is  reversed  at  a  caustic  in  the  wave-number  domain,  where  neigh¬ 
boring  rays  have  the  same  k  but  different  x.  If  there  were  caustics  in  both  domains 
simultaneously,  the  neighboring  rays  would  have  the  same  values  of  x  and  k  and 
would  not  be  distinct. 

The  ray  solution  for  same  quantity  a(xt)  in  the  spatial  domain  has  the  form 

a(x,  t)  =  a0(x)e'[^(x)_®(],  (3) 

with  kfx)  =  Vi/f .  We  assume  for  now  a  simple  time  dependence  e~lajt  with  fixed 
frequency  cu.  The  corresponding  ray  solution  in  the  wave-number  domain  is 

b( k,  t)  =  b0(k)emk)-ro'\  (4) 

with  x(k)  =  — Vk (p.  The  two  phase  functions  (f>  and  ijj  are  related  by  the  Legendre 
transformation  (e.g..  Brown  2000).  Of  more  concern  is  the  relation  between  the 
two  amplitude  functions  ao  and  bo,  which  is  given  by 

ao|Vkx|  =  (5) 

The  Jacobian  Ykx  is  for  the  ray  transformation  xfk)  that  maps  wave  number  to 
position. 

At  a  caustic,  the  ray  transformation  xfk)  is  multivalued.  In  the  case  of  Figure  1, 
two  nonintersecting  neighboring  rays  in  the  wave-number  domain  map  onto  the 
same  spatial  point  x  at  the  caustic,  where  the  Jacobian  vanishes: 

|Vkx|  =  0.  (6) 

Because  x  =  —  Vk</>  (see  Equation  4),  the  Jacobian  can  also  be  expressed  as  (/>kk|, 
a  term  that  appears  in  the  denominator  of  the  amplitude  of  the  stationary-phase 
solution  (e.g.,  Shutts  1998,  equation  65)  and  accounts  for  the  breakdown  of  the 
stationary-phase  method  at  a  caustic.  The  stationary-phase  condition  is  simply  the 
ray  transformation  xfk). 

In  the  following,  we  identify  a[t  and  with  the  wave-action  densities  in  the 
spatial  domain  and  in  the  wave-number  domain,  respectively.  We  continue  to  ignore 
time  dependence  (except  in  the  wave  phase).  We  invoke  the  idea  of  a  narrow 
ray  tube,  which  consists  of  a  set  of  neighboring  rays.  The  wave-action  density 
within  the  ray  tube  is  controlled  by  the  convergence  and  divergence  of  neighboring 
rays  that  make  up  the  ray  tube  (e.g.,  Lighthill  1978).  Conservation  of  wave  action 
implies  that  V  •  (CgCig)  —  0.  Use  of  the  Gauss  divergence  theorem  then  leads  to 
the  relation  for  the  constancy  of  wave-action  flux  through  the  ray  tube  (e.g..  Broad 
1999). 

flo  dx/dt  ■  ndS  —  constant.  (7) 

The  constant  here  (and  in  the  following  three  equations)  is  generally  different  for 
each  ray  tube.  The  group  velocity  is  cg  —  dx/dt ,  and  n  dS  is  a  directed  area 


238  BROUTMAN  ■  ROTTMAN  ■  ECKERMANN 


element  spanning  the  ray  tube  (e.g.,  Lighthill  1978,  figure  89).  An  analogous 
solution  exists  for  the  wave-number  domain: 

b\d\n/dt  •  ndS  =  constant,  (8) 

where  n  dS  is  a  directed  area  element  spanning  the  ray  tube  in  the  wave-number 
domain.  We  return  to  Equation  8  in  Section  5.3. 

Two  cases  of  Equation  7  are  used  commonly.  In  the  first  case,  n  is  parallel  to 
c g,  and  dS  is  the  cross-sectional  area  of  the  ray  tube.  The  left  side  of  Equation  7 
reduces  to  afagldS,  as  Lighthill  (1978,  p.  321)  notes. 

In  the  second  case,  n  is  directed  vertically,  along  the  z  axis.  Then  dS  measures 
the  area  of  a  horizontal  slice  through  the  ray  tube,  and  wave-action  conservation 
becomes 


a^CgiJi  —  constant.  (9) 

The  Jacobian  J\  —  d(x,  y)/d(xo,  yo)  is  taken  at  fixed  z  and  is  proportional  to  dS. 

The  vertical  component  of  the  group  velocity  is  cg 3,  and  the  coordinates  x$,  yo 
refer  to  a  reference  position.  [For  a  horizontally  uniform  background,  it  is  some¬ 
times  convenient  to  use  k,  l  instead  of  xq,  vo,  as  seen  in  Shutts  (1998)  and  Broad 
(1999)]. 

If  J\  is  constant,  we  are  left  with  the  constancy  of  c^a^,  the  vertical  flux  of  wave 
action.  This  has  probably  been  the  most  widely  used  equation  for  wave-amplitude 
calculations  in  internal-wave  models.  The  assumption  (in  addition  to  the  neglect 
of  time  dependence)  is  that  the  horizontal  divergence  of  the  rays  can  be  neglected. 

So  far,  we  have  considered  representations  in  spatial  coordinates  and  in  wave- 
number  coordinates.  It  is  sometimes  convenient  to  mix  the  coordinates,  e.g.,  to 
combine  the  vertical  spatial  coordinate  z  with  the  horizontal  wave-number  coor¬ 
dinates  k,  l.  All  of  the  above  relations  generalize  to  this  case  in  a  natural  way.  For 
instance,  the  expression  analogous  to  Equation  9  is 

qo2cg3  Ji  =  constant.  (10) 

Here  the  Jacobian  J2  —  3 (k,  l)/d(ko ,  Iq)  forreference  values  ko,  lo;  the  wave-action 
density  in  (k,  l,  z)  is  denoted  by  qo1.  This  form  is  especially  useful  for  a  horizontally 
uniform  medium  because  the  horizontal  wave  numbers  are  then  constant  along  the 
ray,  leaving  J2  —  l.  In  Section  4,  we  discuss  the  mapping  of  the  ray  solution 
associated  with  Equation  10  into  a  spatial  solution  by  inverse  Fourier  transform. 

Finally,  we  discuss  ray  formulations  for  the  general  case  that  includes  full  time 
dependence  as  well  as  full  spatial  dependence.  The  general  expression  for  wave- 
action  conservation,  in  the  form  of  the  ray  equation,  is 

dA/dt  =  -AV  •  cg.  (11) 

Here  we  use  the  notation  A  (previously  al  but  now  allowed  to  be  time  dependent) 
for  the  wave-action  density  in  the  spatial  domain.  Chain-rule  differentiation  of 
V  ■  cg  leads  to  an  expression  involving  the  wave-number  gradient  tensor  Vk. 
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A  separate  ray  equation  for  Vk  can  be  derived  (Hayes  1970),  but  elements  of 
Vk  diverge  whenever  the  ray  meets  a  caustic.  So,  near  a  caustic,  Hayes  (1970) 
suggests  reformulating  the  ray  equations  in  terms  of  V^x,  the  inverse  of  Vk.  This 
is  in  effect  a  mapping  to  the  wave-number  domain,  which  avoids  the  caustic  in 
the  spatial  domain.  Generally,  there  are  also  caustics  in  the  wave -number  domain, 
where  V^x  diverges,  so  a  ray-tracing  scheme  of  this  kind  would  need  to  alternate 
between  formulations  based  on  the  two  tensors  Vk  and  Vrx. 

A  single  formulation  is  obtained  from  the  parametric  representation  x(a,  t) 
and  k(a,  t).  Here  a,  which  has  the  same  dimension  as  x,  is  a  label  for  each  ray, 
for  instance  its  initial  position.  The  quantity  of  interest  for  the  calculation  of  the 
spatial  wave-action  density  A  is  the  Jacobian  J  —  |V„x|,  taken  at  fixed  f,  which 
measures  the  changing  volume  of  an  element  of  the  ray  tube  advected  along  the 
ray  at  the  local  group  velocity.  At  a  caustic,  J  vanishes  so  the  ray  integration  can 
proceed  through  the  caustic  without  dealing  with  singularity  quantities  such  as  A. 
The  wave-action  density  is  computed  at  positions  before  or  after  the  caustic  using 
the  constancy  of  AJ  along  the  ray,  though  as  Brown  (2000)  states,  this  does  not 
alter  the  fact  that  the  ray  solution  breaks  down  near  the  caustic.  To  obtain  /,  we 
need  in  the  general  case  a  ray  equation  for  the  nonsymmetric  tensors  Vax  and  Vak 
(Hayes  1970,  White  &  Fornberg  1998): 

c/Vax/c/f  =  Vax  •  £2xk  +  Vak  •  S2kk  (12) 

c/Vak/c/f  =  — Vax  •  —  Vak  ■  f2kx.  (13) 

These  are  derived  from  Equations  1  and  2  by  applying  the  operator  Va.  Note  that 
Va  and  d/dt  commute,  unlike  V  and  d/dt.  Similar  equations  have  been  used  to 
assess  the  stability  of  trajectories  in  a  Hamiltonian  system  (e.g.,  Gutzwiller  1990, 
p.  88).  In  some  applications  it  is  enough  to  know  the  total  amount  of  wave  action 
carried  along  with  a  group  of  waves.  This  amount  is  not  dependent  on  the  focusing 
of  rays  within  the  group  and  is  constant  following  the  group: 

/  A(x,  t)  dx  —  constant.  (14) 

Jvu) 

The  integral  is  taken  over  a  volume  V(t ),  whose  boundaries  move  at  the  local  group 
velocity.  For  an  infinitesimally  sized  V{t),  A  can  be  removed  from  the  integral, 
and  Equation  14  is  equivalent  to  the  constancy  of  AJ  along  the  ray.  An  analogous 
result  holds  for  the  wave-number  domain: 

/  fi(k,  t)dk  =  constant.  (15) 

JK(t) 

Here  B  (previously  but  now  allowed  to  be  time  dependent)  is  the  wave-number 
density  of  wave  action,  and  the  integral  is  taken  over  a  wave-number  volume  IC(t ) 
that  moves  with  the  rays  in  the  wave-number  domain.  For  a  uniform  medium,  1C 
has  fixed  size,  independent  of  time,  because  k  is  constant  along  the  ray.  For  more 
on  volume  integrals  of  this  type,  see  Btihler  et  al.  (1999). 
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In  the  rest  of  the  paper  we  specialize  to  internal  waves.  We  ignore  non-Boussi- 
nesq  effects,  which  are  important  in  atmospheric  applications  but  not  for  our 
discussion  of  ray  formulations.  We  write  the  dispersion  relation  in  the  form 


«  =  &  +  k  U,  (16) 

where  U(x,  t)  is  the  background  flow  and  a>  is  the  intrinsic  frequency: 

cb  =  {k2hN2  +  m2f2)l,2/{k2  +  m2)‘/2  (17) 

The  magnitude  of  the  horizontal  wave  number  is  k/?  =  (k2  +  l2)^2.  The  inertial 
(Coriolis)  frequency  is /and  the  mean  buoyancy  frequency  is  N.  In  the  hydrostatic 
limit  without  Coriolis  effects,  the  dispersion  relation  reduces  to 

a  ±ki,N/m.  (18) 

For  internal  waves  the  wave-energy  density  E  —  Aco  is  related  to  the  vertical 
displacement  amplitude  rj q  by 

E  =  (^2  +  f2m2/k2) ,  (19) 

where  po  is  the  mean  density. 


3.  CAUSTICS 

The  difference  between  ray  theory  and  linear  theory  is  that  the  former  assumes 
slowly  varying  waves.  The  waves  are  not  slowly  varying  in  the  vicinity  of  a  caustic, 
defined  as  the  line  or  surface  containing  points  where  neighboring  rays  intersect 
each  other.  The  simplest  caustic  to  analyze  is  flat,  and  correctable  with  an  Airy 
function.  The  presence  of  caustic  curvature  complicates  the  numerical  implemen¬ 
tation  of  the  caustic  correction.  Note  that  the  Airy  function  solution  for  the  curved 
caustic  pictured  in  Lighthill  (1978,  figure  98)  depends  on  the  caustic  curvature 
through  the  third  derivative  terms  in  Lighthill’s  equation  381.  To  determine  the 
caustic  curvature  we  need  either  to  advect  second  derivatives  of  the  wave  number 
along  the  ray  (in  addition  to  the  first  derivatives  in  Equation  13,  or  to  trace  enough 
rays  to  map  out  the  shape  of  the  caustic.  Neither  appraoch  is  ideal  for  practical 
ray  tracing.  In  addition,  other  types  of  caustics  occur  that  are  not  treatable  with  an 
Airy  function. 

For  example,  consider  the  propagation  of  internal  waves  generated  by  flow  over 
a  mountain.  Such  mountain  waves  (reviewed  by  Wurtele  et  al.  1996  and  Baines 
1995)  often  grow  to  large  amplitudes  by  the  time  they  reach  the  stratosphere.  Their 
dissipation  is  important  for  driving  stratospheric  winds,  so  the  entire  process  of 
mountain-wave  generation,  propagation,  and  dissipation  needs  to  be  parameterized 
in  circulation  models  that  cannot  adequately  resolve  the  mountain  waves.  An  early 
effort  by  Palmer  et  al.  (1986)  was  based  on  a  simple  scheme  designed  by  Lindzen 
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Figure  2  Two  examples  of  an  extreme  breakdown  of  ray  theory.  The  internal  waves 
are  generated  by  (a)  flow  over  mountain  and  (b)  a  localized  source  of  fixed  frequency. 
All  rays  are  confined  to  the  vertical  axis  in  (a)  and  to  the  surface  of  the  cone  in  ( b ). 


(1981).  If  we  try  to  improve  Lindzen’s  scheme  using  a  standard  ray  method,  the 
resulting  ray  prediction  for  the  mountain-wave  amplitudes  is  completely  useless. 

Lindzen  treats  the  mountain  waves  as  a  single  hydrostatic  wave  train  that  prop¬ 
agates  directly  upward  above  each  mountain.  The  ray  solution,  on  the  other  hand, 
includes  rays  for  a  specttum  of  wave  numbers,  but  in  this  case  every  ray  follows 
the  same  path  directly  upward  from  the  center  of  the  mountain  (Figure  2a).  There 
is  no  horizontal  propagation  because  the  intrinsic  horizontal  group  velocity  of 
the  mountain  waves  is  directed  upwind  and  is  perfectly  negated  by  the  horizontal 
background  wind.  The  ray  prediction  for  the  vertical  displacement  amplitude  is 
infinite  on  the  vertical  axis,  where  there  are  infinitely  many  rays  at  each  point,  and 
zero  at  points  off  the  vertical  axis,  where  there  are  no  rays. 

A  similar  problem  affects  Lighthill’s  theory  of  waves  emitted  from  a  localized 
source  (Lighthill  1978,  section  4.9).  The  source  has  fixed  frequency,  and  the  back¬ 
ground  is  stationary  and  uniform.  Lighthill  originally  developed  this  theory  for 
other  types  of  waves.  Its  application  to  internal  waves  breaks  down  because  all 
rays  are  constrained  to  the  surface  of  a  cone  (Figure  2b).  The  ray  prediction  for  the 
vertical  displacement  amplitude  is  infinite  at  points  on  the  cone,  where  there  are 
infinitely  many  rays  at  each  point,  and  zero  at  points  off  the  cone,  where  there  are 
no  rays.  This  cone  is  an  example  of  a  structurally  unstable  caustic,  i.e.,  the  caustic 
can  be  eliminated  by  a  mere  perturbation.  For  example,  moving  the  source  relative 
to  the  background,  even  at  the  slightest  speed,  spreads  the  rays  off  the  cone  and 
makes  the  ray  prediction  finite.  This  is  exactly  what  Lighthill  does  in  section  4.12 
of  his  book.  But  a  finite  ray  solution  is  still  inaccurate  if  the  rays  do  not  separate 
sufficiently,  and  so  the  question  arises:  How  fast  do  we  have  to  move  the  source 
to  get  an  accurate  ray  solution? 

For  the  mountain-wave  problem,  the  vertical  line  of  rays  above  the  mountain 
is  also  a  structurally  unstable  caustic.  Adding  the  slightest  nonhydrostatic  effects 
or  Coriolis  effects  is  enough  to  spread  the  rays  downwind  from  the  vertical  axis 
and  to  give  finite  ray  amplitudes.  But  again,  unless  these  effects  are  significant,  the 
ray  solution  in  the  important  region  directly  over  the  mountain,  at  all  heights,  is 


242 


BROUTMAN  ■  ROTTMAN  ■  ECKERMANN 


erroneously  large  because  the  rays  are  not  well  separated.  Note  that  the  only  two- 
dimensional  ray  solution  for  mountain  waves  given  in  Baines  (1995,  p.  242)  uses 
delta-function  topography,  which  generates  a  strongly  nonhydrostatic  response. 
Similar  problems  occur  in  the  three-dimensional  case.  Although  some  rays  disperse 
downwind  of  the  mountain  in  three  dimensions,  as  part  of  a  ship-wave  pattern  of  lee 
waves,  other  rays  remain  strongly  focused  above  the  mountain.  (See  the  singularity 
in  the  stationary-phase  approximations  of  Smith  1980  and  of  Shutts  1998,  and  see 
figure  lb  of  Broutman  et  al.  2002). 

The  feature  that  distinguishes  the  Airy  function  caustic  from  the  caustics  men¬ 
tioned  above  is  the  number  of  ray  intersections.  In  the  above  cases,  an  infinite 
number  of  neighboring  rays  intersect  at  each  point  on  the  caustic.  For  the  Airy 
function  caustic,  only  two  neighboring  rays  intersect  at  each  point  on  the  caustic. 
The  Airy  function  caustic  is  structurally  stable,  a  result  predicted  by  catastrophe 
theory.  For  more  on  caustics  as  catastrophes,  see  Brown  (2000),  Kravtsov  &  Orlov 
(1999),  and  Nye  (1999).  [The  amphidromic  point  in  oceanographic  tidal  maps 
is  another  example  of  a  structurally  stable  singularity  of  wave  theory — see  Nye’s 
book  (Nye  1999).]  In  the  next  section  we  describe  a  method  that  avoids  the  problem 
of  these  structurally  unstable  caustics  without  having  to  worry  about  the  strength 
of  the  perturbation. 


4.  MASLOV'S  METHOD 


We  note  in  the  description  of  Figure  1  that  one  can  avoid  a  caustic  in  the  spatial 
domain  by  mapping  the  rays  to  the  wave-number  domain.  The  mapping  separates 
the  rays  that  intersect  at  a  spatial  caustic.  Maslov’s  method  uses  the  ray  solution 
from  one  domain  to  correct  the  ray  solution  near  caustics  in  the  other  domain. 
To  review  Maslov’s  method,  we  first  consider  the  following  three  possibilities  for 
transforming  from  the  wave-number  domain  to  the  spatial  domain: 


k  —  domain  solution 


transform  x  —  domain  solution 


linear 

ray 

ray 


IFT 

— >•  linear 

ray 

— >  ray 

IFT 

— >  linear  (approx.) 

The  first  procedure  is  the  usual  Fourier-transform  method,  where  IFT  stands 
for  the  inverse  Fourier  transform.  Because  the  ray  solution  is  not  used,  there  are  no 
problems  with  caustics.  But  only  a  restrictive  range  of  applications  can  be  treated 
in  this  way  due  to  the  difficulty  of  finding  the  linear  solution  in  the  wave-number 
domain.  The  second  procedure  starts  with  the  ray  solution  in  the  wave-number 
domain  and  maps  it  to  the  spatial  domain  using  the  ray  mapping  (or  stationary- 
phase  condition)  x(k,  1).  The  ray  mapping  is  multivalued  at  a  caustic,  where  the  ray 
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prediction  gives  infinite  wave  amplitudes.  The  third  procedure  is  a  combination  of 
the  first  two.  It  starts  with  the  ray  solution  in  the  wave-number  domain  and  maps  it 
to  the  spatial  domain  by  IFT.  For  the  moment,  we  assume  that  the  ray  solution  in 
the  wave-number  domain  is  an  accurate  representation  of  the  linear  solution  in  the 
wave-number  domain  (which  rules  out  caustics  in  the  wave-number  domain).  The 
third  procedure  then  gives  an  approximation  to  the  linear  spatial  solution  that  is 
valid  at  all  types  of  caustics  in  the  spatial  domain  without  further  correction.  The 
third  procedure  works  at  spatial  caustics  for  the  same  reason  that  the  first  procedure 
works:  the  IFT  superimposes  all  Fourier  components  to  account  automatically  for 
diffraction  as  needed  near  any  type  of  caustic.  Away  from  the  caustic  there  is  the 
proper  transition  to  the  spatial  ray  solution,  which  appears  automatically  as  the 
stationary-phase  limit  of  the  IFT. 

Note  that  the  IFT  and  the  ray  approximation  do  not  commute.  Taking  the  ray 
approximation  after  the  IFT  (a  further  step  in  the  first  procedure)  yields  the  spatial 
ray  solution  that  breaks  down  at  the  spatial  caustics.  Taking  the  ray  approximation 
before  the  IFT  (the  third  procedure)  yields  a  solution  that  is  valid  at  the  spatial 
caustics,  and  elsewhere,  again  assuming  that  we  start  with  an  accurate  ray  approx¬ 
imation  in  the  wave-number  domain. 

The  third  procedure  is  a  simple  example  of  Maslov’s  method.  Consider  the  case 
of  stationary  mountain  waves  in  a  height-dependent  background.  It  is  convenient 
here  to  use  the  mixed  formulation  k,l,z,  as  in  Equation  10.  Maslov’s  solution  for, 
say,  the  vertical  displacement  i](\)  is  then 

r](x)  =  If  ^rl0(kj,z)ei^m^'>d^eikx+lydkdl.  (20) 


The  term  in  square  brackets  is  the  ray  solution  in  k,  l,  z.  coordinates.  Note  that 
it  has  the  ray-solution  property  that  differentiation  of  the  phase  with  respect  to 
the  independent  variables  k,  l,  z  gives  the  conjugate  variables  —  x,  —y,  m.  For 
differentiation  with  respect  to  z,  this  is  obvious.  For  differentiation  with  respect 
to  k,  /,  the  result  follows  from  =  —dx/dz,  m/  =  —dy/dz  (see  Hayes  1970, 
equation  27a). 

The  amplitude  rj o  in  Equation  20  is  determined  from  conservation  of  wave  action 
in  the  form  of  Equation  10,  and  from  the  lower  boundary  condition  at  z  =  0.  For 
hydrostatic  mountain  waves,  without  the  effects  of  Earth’s  rotation.  Equation  20 
becomes 


rj(x)  = 


h(k ,  I)[m(k,  l,  z)/mQ(k ,  /)]1/2e‘7o  '"(W.zWj  e^+h  dk  dL 


(21) 


Here  h(k,  l )  is  the  Fourier  transform  of  the  mountain,  and  mo  is  the  vertical  wave 
number  at  the  ground  z  —  0.  For  height-dependent  N (z  ),  the  integral  in  Equation  2 1 
should  be  multiplied  by  N(0)/N(z).  Miles  (1969,  equation  4.13),  Shutts  (1998, 
equation  53),  and  Broutman  et  al.  (2002,  equation  27)  derived  solutions  of  this 
type.  None  of  these  studies  accommodates  trapped  waves.  Broutman  et  al.  (2003) 
gives  modifications  for  trapped  waves  and  associated  caustics. 
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Maslov’s  method  also  provides  a  solution  for  internal  waves  radiated  from 
an  oscillatory  source,  the  problem  described  in  the  previous  section.  Lighthill 
(1978,  section  4.10)  treats  this  problem  with  a  Fourier  integral,  but  he  reduces  it 
to  one  dimension  by  applying  the  stationary-phase  method  in  the  other  dimension 
(see  equation  350  in  Lighthill  1978).  That  makes  Lighthill’s  solution  a  far-held 
approximation.  Maslov’s  result  is  valid  in  the  near  held  as  well  as  in  the  far 
held.  There  is  no  distinction  between  the  near  held  and  the  far  held  for  the  k,  l,  z. 
formulation  because  the  rays  are  everywhere  equally  well  separated  by  their  values 
of  k,l.  Lighthill’s  solution  is  also  restricted  to  a  uniform  background  at  rest  with 
respect  to  the  source,  whereas  Maslov’s  solution  applies  to  a  sufficiently  smooth 
but  otherwise  arbitrary  height-dependent  background. 

The  difficult  case  for  Maslov’s  method  is  when  there  are  caustics  in  the  wave- 
number  domain  as  well  as  in  the  spatial  domain.  In  some  cases,  it  is  straightforward 
to  correct  the  caustics  directly  in  the  wave-number  domain,  before  taking  the  IFT, 
as  in  Broutman  et  al.  (2003).  Alternatively,  one  can  apply  the  IFT  to  ray  solutions 
obtained  in  local  regions  of  the  wave-number  domain  where  there  are  no  caustics, 
as  Maslov  showed  with  an  asymptotic  theory.  The  result  of  the  IFT  then  replaces 
the  spatial  ray  solution,  but  only  in  regions  surrounding  the  spatial  caustics.  In 
other  regions,  the  spatial  ray  solution  is  retained.  It  is  not  clear  how  practical  such 
a  procedure  would  be  for  internal  waves,  though  Brown  (2000)  applied  it  success¬ 
fully  to  surface  gravity  waves.  Ziolkowski  &  Deschamps  (1984)  and  Thomson  & 
Chapman  (1985)  discuss  other  applications  of  Maslov’s  method.  Maslov’s  original 
work,  from  the  1960s,  is  summarized  by  Maslov  &  Fedoriuk  (1981). 

5.  APPLICATIONS 

We  now  discuss  a  selection  of  applications,  and  we  continue  to  stress  ray  formu¬ 
lations  rather  than  ray  results.  We  base  the  discussion  on  wave  action,  although 
in  some  models  the  related  quantity  known  as  pseudomomentum  is  of  more  in¬ 
terest  (e.g.,  Warner  &  McIntyre  1996).  The  wave-action  density  is  denoted  by 
A  for  the  spatial  formulation,  and  by  B  for  the  spectral  or  mixed  spatial/spectral 
formulation. 

5.1.  Shear-Generated  Internal  Waves  that  Reach  the  M  esosphere 

Shear  instability  on  the  upper  edge  of  the  jet  stream  leads  to  mixing  patches  whose 
collapse  excites  internal  waves.  Those  internal  waves  that  reach  the  mesosphere 
(at  altitudes  of  50-90  km)  are  potentially  important  for  driving  mesopheric  winds, 
as  discussed  by  Btihler  et  al.  (1999)  and  Biihler  &  McIntyre.  They  used  a  Fourier 
integral  representation  for  the  waves  in  the  near  field  surrounding  the  mixing 
patch,  assuming  a  uniform  background,  and  a  ray  representation  for  the  waves  in 
the  far  field.  The  ray  representation  accounts  for  wave  propagation  through  height- 
dependent  winds.  The  main  concern  is  with  the  total  amount  of  wave  action  that 
reaches  the  mesosphere. 
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Ray  paths  for  this  model  are  plotted  in  Figure  3  using  the  wind  profile  of  Btihler 
&  McIntyre  (1999),  as  indicated  in  the  figure.  The  only  rays  that  can  reach  the 
mesosphere  through  this  wind  are  those  that  leave  the  source  with  an  intrinsic 
group  velocity  in  the  positive  x  and  positive  z  directions.  The  usual  notation  for 
internal  waves  (e.g..  Gill  1982)  is  such  that  these  rays  have  k  >  0  and  m  <  0.  The 
other  rays  are  absorbed  by  critical  layers  or  reflected  by  turning  points  at  altitudes 
below  the  mesosphere. 

The  total  amount  of  wave  action  emitted  by  the  mixing  patch  and  associated 
with  rays  that  have  a  chance  of  propagating  into  the  mesosphere  is  (compare  with 
equation  27  in  Btihler  et  al.  1999) 

+O0  n 

/  B(k)  dk  dl  dm.  (22) 

OO  J  m<0 

Here  B  is  the  wave-action  density  in  the  wave-number  domain.  For  a  uniform 
background,  B  is  independent  of  time,  so  B  and  P  can  be  calculated  from  the 
initial  conditions,  i.e.,  from  the  Fourier  transform  of  the  initial  configuration 
of  the  mixing  patch.  We  have  been  stressing  the  difference  between  the  wave- 
number  and  spatial  formulations  of  ray  theory.  The  spatial  distribution  of  wave 
action  is  highly  time  dependent.  The  spatial  rays  are  given  by  x  =  c gt,  so  the 
spatial  wave-action  density  is  initially  concentrated  at  a  single  point  at  the  cen¬ 
ter  of  the  mixing  patch,  before  dispersing  rapidly  in  all  directions.  The  spatial 
ray  solution  is  not  valid  near  the  mixing  patch  because  the  rays  are  not  sep¬ 
arated  sufficiently.  But  in  wave-number  space,  the  rays  are  separated  by  their 
k  values.  When  Biihler  et  al.  (1999)  calculated  the  near-field  solution  in  Fourier 
space,  they  calculated  the  equivalent  of  the  ray  solution  in  the  wave-number 
domain. 

The  total  wave-action  P  emitted  by  the  mixing  patch  would  be  the  total  wave 
action  received  by  the  mesosphere,  except  that  some  waves  are  reflected  from 
turning  points  before  reaching  the  mesosphere,  and  all  of  the  waves  experience 
viscous  and  radiative  damping,  which  is  important  at  these  altitudes.  These  effects 
need  to  be  taken  into  account,  and  this  is  where  ray  tracing  is  useful. 

Suppose  we  divide  the  spectrum  k  >  0  and  m  <  0  into  contiguous  wave- 
number  sections,  labeled  /C,  for  i  —  1,2,3. ...  If  the  wave-number  sections  are 
small  enough,  k  is  approximately  uniform  within  each  wave-number  section.  We 
can  then  associate  one  ray  and  one  value  of  the  spectral  wave-action  density  B 
with  each  wave-number  section. 

The  wave  action  integrated  over  each  wave-number  section  /C,  is  constant,  as 
expressed  by  Equation  15,  apart  from  damping  effects  that  can  be  modeled  by 


dPi/dt  —  —Pj  ■  [damping  terms],  (23) 

Here  P,  =  lic  it)  If  the  damping  terms  are  parameterized  as  a  function  of 
k,  the  ray  integration  of  Equation  23  is  simple.  The  convergence  of  rays  is  not  a 
concern,  as  it  would  be  in  the  calculation  of  B  itself,  and  caustics  are  irrelevant.  At 
a  caustic  B  diverges  but  the  size  of  the  corresponding  volume  element  /C,  vanishes. 
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The  integral  P,  remains  finite  and  gives  the  proper  ray  prediction  for  the  amount 
of  wave  action  transported  through  the  caustic. 

The  rays  are  integrated  from  the  mixing  patch  to  the  mesosphere,  and  the  total 
amount  of  damping  acting  on  each  P,  is  calculated.  The  damping  factor,  i.e.,  the 
ratio  of  the  final  P,  to  the  initial  P,,  is  then  used  to  multiply  B  in  the  corresponding 
wave-number  element  of  a  discretized  Equation  22.  In  this  way,  the  near-field 
integral  Equation  22  is  modified  to  give  the  far-field  wave  action  received  by  the 
mesosphere. 

The  only  difference  between  the  above  approach  and  that  of  Biihler  &  Mc¬ 
Intyre  (1999)  is  that  the  latter  used  the  phase-space  representation  of  wave  action 
7V(k,  x,  t)  in  place  of  our  wave-number  integral  P, .  Both  JV  and  P,  share  the  prop¬ 
erty  that  for  nondissipative  propagation  they  are  constant  along  the  ray,  unaffected 
by  the  convergence  of  neighboring  rays. 

To  assess  the  validity  of  the  slowly  varying  approximation,  Biihler  &  McIntyre 
(1999)  monitored  the  quantity  m~2dm/dz,  where  dm/dz  —  ( dm/dt)/(dz/dt ). 
This  is  the  best  that  can  be  done  for  the  variables  integrated  in  the  model.  The 
quantity  evidently  approximates  m~2dm/dz,  the  fractional  change  in  m  over  a 
distance  of  m~l  (Lighthill  1978,  equation  139).  The  smallness  of  this  fractional 
change  is  the  appropriate  condition  for  slow  variation  in  certain  one-dimensional 
models.  In  more  than  one  dimension  the  validity  conditions  presumably  involve 
derivatives  of  the  other  wave-number  components,  but  the  general  form  for  the 
conditions  is  not  clear.  Various  conditions  that  we  have  tested,  though  sometimes 
helpful,  do  not  generally  give  a  reliable  indication  of  where  ray  theory  breaks 
down. 

5.2.  Mountain  Waves 

The  simplest  model  of  mountain  waves  is  hydrostatic  and  two  dimensional,  and  it 
results  in  the  worst  possible  breakdown  of  the  slowly  varying  approximation  (see 
Section  3).  All  ray  paths  coincide  on  the  vertical  axis  directly  over  the  mountain. 
An  alternative  is  to  attempt  to  represent  the  average  conditions  over  the  mountain 
with  a  single  ray  tube.  In  the  simplest  arrangement,  the  ray  tube  has  constant  width 
and  is  directed  vertically,  and  the  wave-action  flux  cgjA  is  constant  along  the  ray 
tube,  until  dissipation.  This  idea  has  been  used  in  schemes  for  the  parameterization 
of  mountain-wave  drag  (see  the  review  by  Kim  et  al.  2003)  and  for  operational 
mountain-wave  forecasting  (Bacmeister  et  al.  1994). 

An  improvement  to  this  appoach  is  to  use  several  rays  and  allow  them  to  prop¬ 
agate  laterally  away  from  the  mountain.  The  position  of  each  ray  is  determined 
by  ray  tracing,  but  the  calculation  is  kept  simple  by  preserving  the  vertical  flux  of 
wave  action  cgj  A  for  each  ray,  in  the  absence  of  dissipation.  Some  nonhydrostatic 
and/or  three-dimensional  effects  can  be  incorporated  in  this  way.  Schoeberl  (1985) 
and  Dunkerton  (1981)  give  examples  of  this.  Eckermann  &  Preusse  (1999)  also 
used  this  approach  to  improve  the  forecast  model  of  Bacmeister  et  al.  (1994)  with 
the  ray-tracing  code  of  Eckermann  &  Marks  (1997). 
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The  above  studies  do  not  include  the  effects  of  the  horizontal  divergence  of 
the  rays  on  the  wave  amplitudes.  For  an  infinitesimal  ray  tube,  the  horizontal 
divergence  is  determined  by  a  Jacobian,  as  in  Equation  9.  Shutts  (1998)  and  Broad 
(1999)  made  calculations  using  Equation  9  to  examine  the  approach  of  hydrostatic 
mountain  waves  to  critical  layers  in  three  dimensions.  They  did  not  evaluate  the 
constant  in  Equation  9  but  still  predicted  relative  changes  in  A  along  the  ray. 

An  alternative  is  to  formulate  the  ray  solution  ink,  /  ,  z  and  use  Maslov’s  method. 
Figure  4  shows  an  example  of  Maslov’s  spatial  solution  for  mountain  waves  over 
Scandinavia.  The  calculation  was  made  by  the  authors  in  a  NASA  measurement 
program  during  January  2003.  Vertical  profiles  (assumed  horizontally  uniform) 
for  the  mean  winds  and  the  mean  density  were  obtained  from  a  weather  forecast 
model.  Broutman  et  al.  (2002)  gives  more  details  on  the  calculation  of  Maslov’s 
solution. 

For  the  nonhydrostatic  case,  the  presence  of  trapped  mountain  waves  compli¬ 
cates  the  calculation  of  the  ray  solution  in  both  x,  y,  z  and  k,  l,  z  because  there 
are  turning  points  where  u>  =  N  and  where  cg 3  =  0.  We  are  used  to  thinking  of 
such  turning  points  as  caustics  (e.g.,  Lighthill  1978,  p.  396),  but  this  is  only  true  of 
the  k ,  l,  z  formulation.  The  turning  point  is  not  a  caustic  in  the  spatial  formulation 
because  there  cannot  be  simultaneous  caustics  in  two  different  projections  of  the 
phase-space  rays  (Section  2). 

To  illustrate  this  point.  Figure  5  shows  the  spatial  rays  for  the  same  prob¬ 
lem  presented  in  Figure  1  of  Wurtele  et  al.  (1996),  which  is  also  described  in 
Wurtele  et  al.  (1987)  and  in  Broutman  et  al.  (2003).  The  mountain  is  centered  at 
the  origin,  and  the  wind  flows  from  left  to  right,  increasing  linearly  with  height. 
The  rays  do  not  intersect  at  a  turning  point,  and  hence  the  turning  point  is  not  a 
caustic  in  the  x,y,z  formulation.  The  rays  encounter  caustics,  which  appear  as 
approximately  straight  lines  that  slope  upwards  from  the  origin.  Note  that  each  ray 
reflects  from  its  turning  point  at  a  position  that  is  slightly  to  the  right  of  the  nearest 
caustic.  For  more  on  the  caustics  in  this  particular  problem,  see  Broutman  et  al. 
(2003). 

5.3.  M  odels  of  I  nternal- Wave  Spectra 

We  now  consider  models  that  combine  ray  methods  with  a  statistical  representation 
of  the  wave  field.  We  start  with  a  case  from  the  ocean:  the  refraction  of  short  internal 
waves  by  a  spectrum  of  longer  internal  waves.  Using  ray  methods  for  this  problem 
began  in  earnest  with  the  preliminary  study  of  Henyey  &  Pomphrey  (1983),  and 
continued  with  Flatte  et  al.  (1985),  Henyey  et  al.  (1986),  and  Sun  &  Kunze  (1999). 
(See  also  the  review  by  Muller  et  al.  1986).  These  studies  implement  Monte  Carlo 
ray  tracings  involving  the  Garrett-Munk  model  spectrum  (Garrett  &  Munk  1979), 
which  approximates  measurements  from  the  ocean  and  which  is  used  in  two  ways: 
to  set  the  amplitudes  of  the  background  long  waves,  and  to  set  the  initial  conditions 
for  the  short  waves.  The  idea  is  to  duplicate  the  Garrett-Munk  model  in  the  initial 
conditions,  for  short  waves  of  relatively  large  wavelengths,  and  then  to  see  if 
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the  ray  tracing  duplicates  the  Garrett-Munk  model  at  smaller  scales,  down  to 
dissipative  wavelengths.  Also  of  interest  is  the  flux  of  wave  action  through  the 
spectrum,  which  is  used  to  predict  ocean-mixing  rates  resulting  from  internal-wave 
dissipation. 

The  long-wave  background  is  variable  in  all  three  spatial  dimensions  and  in 
time.  This  is  not  a  serious  complication  for  calculating  ray  paths,  but  it  is  a  serious 
complication  for  calculating  the  wave-number  density  or  spatial  density  of  wave 
action  along  the  ray.  The  ray  tracing  would  require  the  initialization  and  integration 
of  the  full  set  of  ray  equations  in  Equations  12  and  13,  and  would  undoubtedly 
lead  to  frequent  occurrences  of  caustics. 

Flatte  et  al.  (1985),  Henyey  et  al.  (1986),  and  Sun  &  Kunze  (1999)  simplified 
the  wave-amplitude  calculation  and  eliminated  caustics  by  defining  the  ray  tube 
statistically.  They  assumed  that  the  statistics  represented  an  internal-wave  spectrum 
that  was  stationary  and  horizontally  isotropic. 

To  see  how  this  works,  consider  the  form  of  wave-action  conservation 

B(k/,,  m)dki,/dtAm  =  constant,  (24) 

where  B  is  the  spectral  wave-action  density.  This  is  a  special  case  of  Equation  8 
for  the  constancy  of  the  wave-action  flux  b^dk/dt  hdS  along  a  ray  tube  in  the 
wave-number  domain.  To  obtain  Equation  24  from  Equation  8,  n  is  the  direction 
of  the  horizontal  wave-number  axis,  and  the  width  of  the  ray  tube  dS  is  the  vertical 
wave-number  variation  Am  across  the  ray  tube.  It  is  then  assumed  that  each  term 
in  Equation  24  can  be  represented  by  its  averaged  value. 

For  example.  Sun  &  Kunze  (1999)  used  this  approach  to  estimate  the  flux  of 
wave  energy  to  short  dissipative  scales.  Initial  conditions  are  specified  at  the  rela¬ 
tively  large  scale  for  the  short  waves  of  1-km  horizontal  wavelength:  the  Garrett- 
Munk  spectrum  sets  the  initial  average  for  B ,  the  discretization  of  the  Garrett-Munk 
spectrum  determines  the  initial  Am,  and  various  estimates  are  used  to  set  the  initial 
average  value  for  dk/,/dt  (see  Sun  &  Kunze,  p.  2912).  This  determines  the  con¬ 
stant  on  the  right  side  of  Equation  24,  i.e.,  the  wave-action  flux  for  each  ray  tube. 
The  wave-energy  flux  is  then  calculated  at  a  small  dissipative  scale  by  multiplying 
the  wave-action  flux  by  the  short-wave  intrinsic  frequency  a>  at  that  small  scale. 
Here  a>  is  obtained  by  tracing  individual  rays  through  realizations  of  the  Garrett- 
Munk  background.  The  small  dissipative  scale  was  chosen  to  be  5-m  vertical 
wavelength. 

Figure  6  illustrates  another  point  about  these  ocean  ray  models.  The  ray  paths 
shown  in  this  figure  follow  short-wave  propagation  through  a  vertically  localized 
packet  of  inertia  waves  centered  in  the  middle  of  the  plot.  A  numerical  solution  for 
the  short  waves  is  also  shown.  When  this  figure  was  first  published  in  Broutman 
et  al.  (1997),  the  main  interest  was  in  the  initial  encounter  of  the  short  waves 
with  the  inertia-wave  packet,  at  times  just  after  one  inertia  period.  The  permanent 
upturn  of  the  rays  after  one  pass  through  the  inertia-wave  packet  and  the  absense 
of  critical  layers  for  the  short  waves  were  noted  as  special  features  of  refraction 
by  time-dependent  shear. 


depth  in  inertia  wavelengths 
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Figure  6  Ray  paths  and  a  numerical  solution  for  short  waves  propagating  through  an 
inertia-wave  packet.  From  Broutman  et  al.  (1997). 


Here  we  shift  attention  to  earlier  times.  The  short  waves  wrap  around  the  pe¬ 
riodic  computational  domain  and  repeatedly  encounter  the  inertia-wave  packet. 
The  rays  show  signs  of  chaotic  behavior,  as  they  are  likely  to  do  in  the  models 
of  Flatte  et  al.  (1985),  Henyey  et  al.  (1986),  and  Sun  &  Kunze  (1999).  Henyey 
et  al.  (1986,  section  2)  noted  that  the  individual  rays  in  their  calculations  are  very 
chaotic.  However,  the  statistics  that  they  derive  from  the  individual  rays  appear 
to  be  very  stable,  and  in  good  agreement  with  measurements.  It  may  seem  para¬ 
doxical  that  chaotic  rays,  with  the  sensitivities  in  tracing  them,  can  yield  stable 
and  reliable  results  for  the  short-wave  statistics,  but  this  has  been  shown  to  occur 
in  other  studies  of  ray  chaos,  for  example  in  acoustics  and  semiclassical  physics 
(e.g..  Brown  et  al.  2003). 

Atmospheric  refraction  models  have  been  developed  with  similar  aims  of  pre¬ 
dicting  internal-wave  spectra  and  dissipation  rates.  Warner  &  McIntyre’s  (1996, 
1999,  2001)  models  are  intended  to  be  used  for  wave -drag  parameterization  in 
general  circulation  models  and  are  thus  constrained  by  computational  costs  to  a 
much  simpler  design  than  the  ocean  models  described  above.  For  example,  the  re¬ 
fraction  of  short  waves  by  long  waves  is  ignored,  and  only  the  vertical  variability 
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in  the  background  winds  is  taken  into  account  for  the  refraction.  The  rays  are  thus 
not  chaotic,  but  the  wave  field  is  still  represented  statistically  with  a  generalized 
wave  packet,  defined  as  a  discrete  element  of  the  spectrum.  Note  also  the  use  of 
wave-number  (and  frequency)  coordinates  in  their  ray  formulation. 

Hertzog  et  al.  (2002)  and  Souprayen  et  al.  (2001)  give  another  statistical  wave 
formulation.  They  use  a  phase-space  representation  to  study  the  refraction  of  short 
internal  waves  by  a  spectrum  of  long  internal  waves  in  the  atmosphere.  They 
include  the  full  space  and  time-dependence  of  the  long  waves  in  the  ray  tracing. 
The  rays  have  numerous  caustics  when  expressed  in  the  spatial  domain  (appendix 
C  of  Hertzog  et  al.  2002),  but  as  noted  earlier  the  phase-space  formulation  is 
free  of  caustics  and  the  phase-space  density  of  wave-action  Af( k,  x,  t)  is  constant 
along  the  ray.  Estimates  of  the  energy  spectrum  in,  say,  ( m ,  z),  are  obtained  by  an 
integration  of  o>N  over  k,  l,  x,  y  (see  equation  6  of  Hertzog  et  al.  2002). 

In  models  such  as  these,  the  rays  wander  quasi-randomly  through  large  por¬ 
tions  of  the  allowable  phase  space.  The  projection  of  the  rays  onto  the  spatial, 
wave-number,  or  mixed  domains  results  in  a  dense  concentration  of  caustics,  and 
Maslov’s  method  is  not  practical.  The  best  hope  of  dealing  with  caustics  in  this  sit¬ 
uation  is  to  smooth  over  them.  The  integration  preformed  by  Hertzog  et  al.  (2002) 
smoothes  the  caustics  and  seems  to  be  similar  in  some  respects  to  the  treatment 
described  in  Berry  (1983,  section  3.4).  Berry  also  gives  a  useful  discussion  of  the 
representation  of  a  wave  field  in  phase  space. 

5.4.  Other  Applications 

When  the  background  is  time  varying,  the  ray-tube  equations  (Equations  7-10) 
are  not  applicable.  Some  models  have  been  developed  for  a  time-varying  but 
spatially  uniform  background,  e.g.,  Lott  &  Teitelbaum’s  (1993)  mountain-wave 
study.  They  computed  ray  and  caustic  solutions  from  the  stationary  phase  and 
Airy  function  limits  of  an  integral  representation.  Ray-tracing  models  that  de¬ 
scribe  time-dependent  short-wave  refraction  by  a  single  long-wave  packet,  or  by 
a  few  long-wave  packets,  include  Sonmor  &  Klaassen  (2000),  Eckermann  (1997), 
Walterschied  (2000),  Zhong  et  al.  (1996),  Thorpe  (1989),  and  Broutman  &  Young 
(1986).  Sonmor  &  Klaassen  (2000)  gave  a  detailed  analysis  of  short-wave  caustics 
resulting  from  long-wave  shear.  Broutman  &  Young  (1986)  used  a  simple  ray  for¬ 
mulation  involving  the  volume  integral  of  wave  action  to  identify  a  mechanism  for 
the  nondissipative  damping  of  the  long  waves  by  the  short  waves.  A  more  detailed 
calculation  appeared  in  Broutman  &  Grimshaw  (1988). 

Horizontally  varying  backgrounds  have  been  treated  with  ray  methods  for  ap¬ 
plications  such  as  internal-wave  propagation  near  fronts  and  vortices  (e.g.,  Kunze 
1985,  Dunkerton  1984,  Hertzog  et  al.  2001).  Often,  important  insights  are  gained 
from  an  inspection  of  the  ray  paths  without  wave-amplitude  calculations  that  would 
be  complicated  by  caustics.  Pringle  &  Brink  ( 1 999)  provide  a  model  with  a  tractable 
ray  and  Airy-function  caustic  calculation  for  internal  waves  over  a  sloping  bottom 
in  the  presence  of  a  horizontally  sheared,  depth-independent  mean  flow. 


RAY  METHODS  FOR  INTERNAL  WAVES 


251 


6.  CONCLU  SION 

We  discussed  ray  models  for  wave-amplitude  calculations  of  internal  waves,  stress¬ 
ing  ray  formulations  rather  than  ray  results  and  practical  implementations  rather 
than  formal  theory.  We  gave  examples  that  use  spatial  coordinates,  wave-number 
coordinates,  and  phase-space  coordinates,  and  that  are  expressed  in  terms  of  the 
local  density  or  the  volume  integral  of  wave  action.  In  most  of  these  cases,  the  ray 
calculation  is  deterministic.  In  some  cases,  the  wave  amplitudes  are  initialized  sta¬ 
tistically  and  then  followed  along  deterministic  ray  tubes  (e.g.,  Warner  &  McIntyre 
1996),  and  in  other  cases  the  wave  amplitudes  and  the  ray  tubes  are  represented 
statistically  (  Hcuyey  et  al.  1986,  Sun  &  Kunze  1999,  Flatte  et  al.  1985). 

The  choice  of  ray  formulation  affects  not  only  the  difficulty  of  the  ray  calcu¬ 
lation,  but  also  the  extent  to  which  the  waves  satisfy  slowly  varying  assumptions. 
For  example,  a  spatial  caustic  can  be  mapped  away  by  changing  some  or  all  of 
the  spatial  coordinates  to  wave-number  coordinates.  One  formulation  rarely  suits 
the  entire  problem:  the  initialization,  the  ray  tracing,  the  correction  of  caustics, 
the  application  of  dissipative  schemes,  and  the  prediction  of  variables  of  interest. 
An  aim  of  this  paper  has  been  to  discuss  how  a  combination  of  formulations  and 
assumptions  has  contributed  to  the  development  of  practical  ray-tracing  schemes. 
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Figure  3  Ray  paths  for  the  Biihler-Mclntyre  model  (Section  5. 1)  of  internal-wave  prop¬ 
agation  from  the  upper  edge  of  the  jet  stream  through  the  stratosphere  and  into  the 
mesosphere. 
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Equation  20  is  evaluated  by  a  trapezoidal  rule. 
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Figure  5  Ray  paths  for  nonhydrostatic  mountain  waves.  The  calculation  corresponds  to 
the  model  in  figure  1  of  Wurtele  et  al.  (1996).  The  mountain  is  centered  at  the  origin,  and 
the  wind  is  in  the  positive  x  direction  and  increases  with  height. 


